#coding=utf-8

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

nstep = 300
nwalker = 1000
x = [0 for i in range(nwalker)]

def A(xi, xj):
    if abs(xj) < abs(xi): return 1
    else: return np.exp( (xi*xi - xj*xj)/2 )

def onestep():
    for i in range(nwalker):
        step = 0.2*np.random.randint(0,2)-0.1
        if( np.random.random() < A(x[i], x[i]+step) ):
            x[i] += step

def plotNormal():
    x = np.arange(-3,3,0.1)
    y = 1 / np.sqrt(2*np.pi) * np.exp( -x*x/2 )
    plt.plot(x,y, color='red')

fig, ax = plt.subplots()

def update(frame):
    onestep()
    plt.cla(); plt.ylim(0, 1); plt.title("iter = %d"%frame)
    plotNormal()
    print("iter = ", frame)
    plt.hist(x, bins=np.arange(-3, 3, 0.2), rwidth=0.8, density=True )
ani = animation.FuncAnimation(fig, update, range(nstep), interval=100)
#plt.show()
ani.save("Metropolis.gif",fps=60)